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A class of algorithms in discrete space and continuous time for Brownian first passage time 
estimation is considered. A simple algorithm is derived that yields exact mean first passage times 
(MFPT) for linear potentials in one dimension, regardless of the lattice spacing. When applied 
to nonlinear potentials and/or higher spatial dimensions, numerical evidence suggests that this 
algorithm yields MFPT estimates that either outperform or rival Langevin-based (discrete time, 
continuous space) estimates. 



Brownian dynamics is one of the most widespread 
models of temporal evolution for systems displaying 
stochastic behavior 11 . Its popularity stems no doubt 
in part from its simplicity, which allows one to carry out 
analytical work to great lengths, but also from its gener- 
ality, as many dynamical systems ranging from biological 
molecules p] to financial markets [3| are often well ap- 
proximated by this model. 

In general, however, the solution of most Brownian 
problems is not known in closed analytical form, requir- 
ing one to resort to numerical simulations. Almost in- 
variably, such solutions are obtained by discretizing the 
Langevin equation in time, and iterating the ensuing dif- 
ference equation (see e.g. [1]). Here I propose to dis- 
cretize space instead, leaving time continuous. There are 
many ways of going about this procedure, and different 
algorithms can be obtained depending on the desired con- 
text. In this paper I will focus on the design of algorithms 
suited for the computation of mean first passage times 
(MFPT) to a given boundary [5] , which plays a particu- 
larly important role in theories of chemical kinetics [HIIT]. 

To introduce the basic idea behind the present algo- 
rithm, let us focus on the simple one-dimensional prob- 
lem shown in Fig. [T] The illustration depicts a typi- 
cal Brownian trajectory in a first-passage problem from 
a; = to a; = 2 A. This problem is characterized by an 
ensemble of continuous trajectories that start from x = Q 
at f = 0, and cross the absorbing boundary x = 2 A only 
once at some time t = t; t is thus the first passage time 
of the trajectory. Our goal is to design algorithms that 
generate discrete trajectories (blue lines in Fig. [T]) that 
"hop" from site to site so that their MFPT to the ab- 
sorbing boundary approximate that of the original, con- 
tinuous problem. 

The outline of the derivation is as follows. First, the 
mean first passage time of the continuous Brownian prob- 
lem will be recast in terms of two quantities defined on an 
arbitrary lattice, namely conditional mean first passage 
times and splitting probabilities (Eq. This inter- 

mediate result will allow us to design lattice algorithms 
that reproduce the MFPT of the original Brownian prob- 
lem by demanding that their conditional MFPTs and 




FIG. 1: A continuous Brownian trajectory (wiggly black line) 
and its discrete counterpart (straight blue lines), illustrating 
a first passage problem from a; = to the absorbing boundary 
aX X = 2A. The discrete states are labeled s = 0, ±1, ±2, and 
etc, corresponding to a: = 0, ±A, ±2A, and etc. For both 
types of trajectories, the total first passage time r is the sum 
of the conditional first passage times r(si+i|si) from state Si 
to the next state Si+i, where "conditional" means that the 
particle did not cross the other adjacent state before crossing 
Si+i [5]. The discrete trajectories are constructed so that 
their MFPT (t) is the same as that of the original Brownian 
trajectories (see text). 



splitting probabilities be equal to those of the Brown- 
ian problem. As these two quantities are generally not 
algebraic for nonlinear potentials, they will be evaluated 
based on a linear approximation in the region delimited 
by the nearest neighbor sites (Eqs. (|5| and ([6])). By addi- 
tionally demanding the sites to be uniformly spaced, this 
will allow us to write a generic rate equation that can be 
generalized to higher dimensions (Eqs. ^ and (lOl). Fi- 
nally, this rate equation is simulated by standard means, 
e.g. using Gillespie's algorithm [5]. 

Going back to Fig.[l] we see that the first passage time 
of any Brownian trajectory can be decomposed as a sum 
of intermediate times T{si+i\si). Thus, the mean first 
passage time from state si to state sm in the restricted 
ensemble of trajectories that pass through a given time- 
ordered sequence of states — {si, S2, . . . , sat} is 
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The quantity (r(si-(_i|si)) is the conditional mean first 
passage time from state Si to state Si+i, where the term 
"conditional" means that the particle is not allowed to 
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pass through the other adjacent state [S] (e.g. t(1|0) is 
the first passage time from s = to s = 1, conditional on 
not passing through s = —1). Note that the individual 
terms of this sum depend only on the present and next 
states, Si and s^+i respectively. This is only true for 
Markovian dynamics, which is assumed to be the case 
for the present Brownian problem. The total MFPT (r) 
is thus obtained by taking the average of (t(s^)) over all 
permissible sequences of states , i.e. 
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where p{s^) is the probability that the particular se- 
quence of states will be realized, and the double sum 
is over all possible sequences of states that take the parti- 
cle from its original position to the absorbing boundary. 

Now, for Markovian dynamics, the probability p{s^) 
can be decomposed as a product of splitting probabilities 
(j), where (j)(si^i\si) is the probability that a particle orig- 
inally at Si will pass through Si_|_i before passing through 
the other adjacent state (e.g. 0(2|1) is the probability 
that the particle originally at s = 1 will pass through 
s = 2 before s — 0). This finally gives the result 
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The main conclusion from this derivation is that the 
MFPT of our Brownian problem is fully specified by the 
splitting probabilities and conditional MFPTs defined on 
an arbitrary lattice (although we have chosen a uniform 
lattice anticipating the development below, this deriva- 
tion is clearly more general). It thus follows that any 
other dynamical system that has the same 4>{s'\s) and 
{t{s'\s)) for all adjacent sites s' ,s as the original Brown- 
ian problem also has the same MFPT (r). In turn, this 
suggests that the design of MFPT algorithms on a lattice 
should focus on reproducing as closely as possible these 
two quantities from the original Brownian problem. 

For one dimensional Brownian problems, both 0(s ± 
l|s) and (r(s ± l|s)) can be reduced to simple quadra- 
ture |5j. An additional simplification occurs when the 
particle is subject to a linear potential and the lattice is 
uniformly spaced, in which case two things happen: first, 
the integrals reduce to algebraic expressions, and second 
the conditional MFPTs (r(s + l\s)) and (t(s - l|s)) be- 
come coincident and equal to the unconditional mean exit 
time, (r(s)). This second observation allows us to write 
down a rate equation governing the dynamics on the lat- 
tice, which can then be generalized to higher dimensions. 

To be specific, consider a particle evolving according 
to the Smoluchowski equation [1] 



dp 
dt 
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and subject to the linear potential U{x) = ax, where for 
simplicity of notation energy is measured in units of fcsT. 
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FIG. 2: Illustration of a two dimensional implementation of 
Eqs. ([9|-(fT0|. The arrows represent the four outgoing rates 
from s, fc(s±x|s) and fc(s±y|s). The boundary sites lie along 
the vertical line on the right and are highlighted in red. The 
particle "dies" whenever it visits one such site. 



For such one dimensional potentials, the (unconditional) 
mean first passage time [3 |9] from s to the adjacent 
positions s ± 1 is, exactly, 
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while the splitting probabilities are 
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(s±l|s) = 
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Given {t{s)) and (/)(s ± l|s), a lattice rate equation can 
be constructed consistent with these quantities. Indeed, 
consider a kinetic scheme for the states s — 1, s, s-l- 1 with 
outgoing rates from s given by fc(s ± l|s). The lifetime 
(t(s)) in the state s is then {t{s))'^ — k{s— l\s) + k{s + 
l|s), while the splitting probabihties are 0(s ± l|s) = 
(T(s))fc(s ± l|s). Solving these equations for the rates 
and using the results for linear potentials above, we get 



k{s±l\s) = ± 
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where these rates are to be used in the rate equation 



dpjs; t) 
dt 



k{s\s + l)p{s + l;t) + k{s\s - l)p{s - 1; t) 
-[k{s^l\s) + k{s+l\sMs;t). (8) 



Equations ([t]) and ([s]) form the foundation of the pro- 
posed algorithm. For linear potentials in one spatial di- 
mension, the algorithm yields exact MFPTs. For non- 
linear potentials, a is to be replaced by the slope of the 
potential at the position corresponding to site s (local 
linear approximation). In higher dimensions, the rate 
equation Eq. ([s]) can be generalized by taking the rates 
along each coordinate to be the one dimensional result 
already derived. Thus, the general form of our rate equa- 
tion takes the form 



dpjs; t) 
dt 



J2 [k{s\s')p{s';t)-k{s'\s)p{s;t)], (9) 
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where the sum is over the nearest neighbors of s, and 
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/c(s±z|s) = ±- 
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In this last equation, z is a unit basis vector along any 
of the coordinates, and Uz{s) is the partial derivative of 
the potential with respect to that coordinate evaluated 
at the position corresponding to the site s. For simplic- 
ity, cartesian coordinates and square lattices are being 
assumed (see Fig.[2|. 

Before discussing how to simulate the above rate equa- 
tion, let us check that in the continuum limit we are ex- 
actly solving the Smoluchowski equation (Eq. Q). When 
A is small in Eq. (10 1, we have to leading order 
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Substituting these rates into Eq. (|9| and mapping finite 
differences into differential operators, we indeed obtain 
Eq. This shows that, although our method was de- 
signed with MFPT estimation in mind, the ensuing al- 
gorithm actually generates exact trajectories in the con- 
tinuum limit, much like the Langevin algorithm becomes 
exact when the time step goes to zero. 

The simulation of Eqs. ^ and ([lO I can be performed 
by means of Gillespie's celebrated algorithm 8 . Accord- 
ing to this method, one starts in a given state s and draws 
an exponentially distributed random number t with mean 
equal to the reciprocal of the sum of the outgoing rates, 
i.e. {t)~^ = J2s'=nn k{s'\s). This is the lifetime of the 
particle in the state s. A decision is then made as to 
which site among the nearest neighboring states the par- 
ticle is going next. This is done by assigning the statisti- 
cal weight w(s') = k{^'\s)/J2s'=nn k{s'\s) to each neigh- 
boring site s', and choosing one such site with probability 
w{s'). The particle then moves to this chosen site, and 
the procedure above is repeated until the particle reaches 
a boundary site. The sum of the times t until this cri- 
terion is satisfied is then the first passage time to the 
boundary. 

The above algorithm has been tested on model prob- 
lems in one and two spatial dimensions (Figures [s] and |4] 
respectively) . For comparison, the overdamped Langevin 
algorithm was also simulated |1] 

y:{t + At) = x(t) - DAtVU{x{t)) + V2DAtg. (12) 

Here At is the time step, and g is a Gaussian random vec- 
tor of zero mean and unit variance. Whenever available, 
numerically exact results are also reported to illustrate 
the correctness of the algorithm in the continuum limit. 
As force computation is the main bottleneck in most sim- 
ulations, the main figure of merit considered was the av- 
erage number of force evaluations per trajectory required 
to achieve a given accuracy level. In all test cases consid- 
ered, the present algorithm requires considerably fewer 
force evaluations than the Langevin algorithm, although 
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FIG. 3: Numerical results in one dimension, comparing 



Langevin's algorithm (Eq. ( 12 1, red squares) with the present 
algorithm (Eqs. ([7|-([8]l, blue circles). Top: Mean first pas- 
sage time from a; = to a; = 1 for the linear potential 
U{x) = —X. The exact result obtained by analytical integra- 
tion [9; is (r) = 1 (dashed line). Bottom: MFPT from x = 
to X — \/6 for the harmonic potential U{x) = a;^/2. The "ex- 
act" result obtained by numeric quadrature [9] is (r) — 24.324 
(dashed line). For both problems D = 1. The error bars are of 
the size of the symbols and hence not shown. For the Langevin 
algorithm, the symbols correspond to decreasing values of the 
time step At, from left to right (e.g. At = 0.5,0.25,0.125, 
etc). For the present algorithm, the symbols correspond to 
increasing number n of lattice points between the starting 
point and the boundary, from left to right (e.g. n = 0, 1, 2, 3, 
etc). The average number of force evaluations per trajectory 
corresponds to the total number of calls to the function Ux (x) 
divided by the total number of trajectories generated (lO''). 



further experimentation is called for in order to make 
more general conclusions. 

The algorithm considered above is only one of various 
strategies that can be used based on Eq. ^ . A straight- 
forward improvement is to evaluate the conditional MF- 
PTs explicitly for piecewise linear potentials, so that the 
resulting algorithm would be exact for such potentials (as 
opposed to being exact for linear potentials only). This 
too would result in algebraic expressions for (r(s'|s)); 
however its generalization to higher dimensions would be 
non-trivial, as in this case no simple rate equation can be 
written (rate constants imply conditional mean lifetimes 
are the same regardless of the outgoing site). Another 
possibility is to calculate splitting probabilities and con- 
ditional MFPTs for surfaces (instead of lattice points) 
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FIG. 4: Numerical results in two dimensions. Legends and 
parameters are the same as in Figure [3] Top: Escape time 
from the square with corners (0, 0), (1, 0), (1, 0), (1, 1), for 
a free particle starting at (0.5,0.5). The numerically ex- 
act result is (r) = 0.0736714 'Q:. Bottom: MFPT to the 
boundary x = 0.5 for the symmetric double well potential 
U{x,y) = [VSix-lf +y^] x [V3{x + if + y^]. The particle 
starts at the left minimum (—1,0); the other minimum is at 
(1,0). 



surrounding the particle. The advantage in this case is 
that both space and time are treated continuously, and 
the algorithm applies to any number of dimensions (see 
e.g. [rUlfTT] for free-particle implementations). The main 
difficulty here, however, is to find reasonable approxima- 
tions to such quantities when J7(x) ^ 0; to this author's 
knowledge, such algorithms are yet to be designed. 



In summary, in this contribution a new class of al- 
gorithms for the estimation of mean first passage times 
in Brownian dynamics was introduced. In contrast to 
traditional discrete-time (Langevin) methods, these al- 
gorithms treat time continuously and space discretely. 
Perhaps their most distinguishing feature is that they 
can yield exact MFPTs regardless of the lattice spacing 
in some particular cases; for example, the algorithm con- 
sidered above yields exact MFPTs for linear potentials in 
one dimension. Numerical results also suggest that the 
algorithm outperforms Langevin-based estimates in two 
dimensions. Its efficiency in higher dimensions and/or 
more complex geometries is currently under investiga- 
tion. 
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